Elasticity-driven collective motion in active solids and active crystals 
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We introduce a simple model of self-propelled agents connected by linear springs, with no explicit 
alignment rules. Below a critical noise level, the agents self-organize into a collectively translating or 
rotating group. We derive analytical stability conditions for the translating state in an elastic sheet 
approximation. We propose an elasticity-based mechanism that drives convergence to collective 
motion by cascading self-propulsion energy towards lower-energy modes. Given its simplicity and 
ubiquity, such mechanism could play a relevant role in various biological and robotic swarms. 



Animal groups that move together, such as bacterial 
colonies, insect swarms, bird flocks, or fish schools pJE], 
are all examples of biological systems displaying Collec- 
tive Motion (CM). In recent years, the dynamics of such 
systems (referred to here generically as swarms) has been 
the subject of intense research [6H9]. A number of theo- 
retical models have been introduced to study swarms and 
to develop control rules that achieve similar coordinated 
collective dynamics in groups of autonomous robots [3|9]. 
Despite this proliferation of CM algorithms, there is still 
no comprehensive understanding of the underlying mech- 
anisms that can lead a group of self-propelled agents to 
self-organize and move in a common direction. 

The current CM paradigm has been strongly influenced 
by the seminal work of Vicsek et al. [10 , which intro- 
duced a minimal model for flocking, the Vicsek model, 
that has become a referent in the field [7H9] . This model 
describes a group of point particles advancing at a fixed 
common speed, only coupled through alignment interac- 
tions that steer them towards the mean heading direc- 
tion of all particles within a given radius [T0HT2] . In this 
framework, a swarm can be viewed as a group of self- 
propelled spins with aligning interactions, described by 
a variation of the XY-model [13] where spins advance in 
their pointing direction rather than remaining affixed to 
a lattice. In the continuous limit, this system becomes 
a fluid of self-propelled spins that follows the hydrody- 
namic theory developed in [T4HT7] . More recently, other 
models that do not rely on explicit alignment interactions 
have been introduced. In [18 , for example, CM is driven 
by escape-pursuit interactions only; in [19], by inelastic 
collisions between isotropic agents; and in [20] and [21], 
by short-range radial forces coupled to each agent's turn- 
ing dynamics. Given that Vicsek-like algorithms rely on 
explicit alignment rules to achieve CM [5j [TDJ [22j [23], 



it was initially surprising that such systems could self- 
organize without them. While it can be argued that all 
these models include at least an implicit alignment inter- 
action, it remains unclear if they are all driven to CM 
by the same underlying mechanisms and to what extent 
agents must exchange orientation information, either ex- 
plicitly or implicitly, to achieve CM. 

In this Letter, we introduce a CM mechanism that is 
based on a very different paradigm: the emergence and 
growth of regions of coherent motion due to standard 
elasticity processes. We explore this mechanism by in- 
troducing a simple two-dimensional Active Elastic Sheet 
(AES) model with spring-like interactions between neigh- 
boring agents and no explicit alignment, which describes 
what we refer to as an active solid or an active crystal. 

We define the AES model as a system of N agents on a 
two-dimensional plane, where the position xi and orien- 
tation 0i of each agent i follow the overdamped equations 
of motion 
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Here, vq is the forward biasing speed that induces 
self-propulsion (injecting energy at the individual par- 
ticle level), hi and hf- are two unit vectors point- 
ing parallel and perpendicular to the heading direction 
of agent z, and parameters a and f3 are the inverse 
translational and rotational damping coefficients, respec- 
tively. The total force over agent i is given by Fi — 
^Zj^Si (~ k/Uj) (\xi — xj I — lij), a sum of linear spring- 
like forces with equilibrium distances 1^ and spring con- 
stants k/kj. Each set Si contains all agents interacting 
with agent i and remains fixed throughout the integra- 
tion. This system is thus akin to a spring-mass model 
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of elastic sheet [24] where masses are replaced by self- 
propelled agents that turn according to Fi • hj- and move 
forward or backwards following Fi • hi and their self- 
propulsion. We include actuation noise (fluctuations of 
the individual motion) by adding Dq to the heading 
angle, where Dq is the noise strength coefficient and 
a random variable with standard, zero-centered normal 
probability distribution of variance 1. We include sens- 
ing noise (errors in the measured forces) by adding D r £ r 
to Fi, with D r the noise strength coefficient and £ r a ran- 
domly oriented unit vector. The degree of alignment in 
the system is monitored through the usual polarization 
order parameter 
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If all agents are perfectly aligned, we have ip = 1. If 
agents are instead randomly oriented or rotating about 
the group's bary center, we have ijj = 0. 

The AES model was designed to study CM under con- 
ditions that are in many ways opposite to those in the 
Vicsek model. While the only information that agents 
exchange there is their relative heading angle, here they 
only sense their relative positions. While changing in- 
teracting neighbors over time has been shown to be nec- 
essary there for achieving long-range order [TTJ [14j [25] , 
here virtual springs connect the same agents through- 
out the dynamics. Furthermore, while both models de- 
scribe overdamped systems, these are implemented dif- 
ferently. In the Vicsek model, particles switch instanta- 
neously to the desired heading angle; in the AES case, we 
integrate instead the overdamped angular equation (J2|, 
which turns out to be necessary here for achieving CM. 

We integrate Eqs. ([!]) and Q numerically using a stan- 
dard Euler method. All simulations below are carried out 
with a 0.01, (3 0.12, v 0.002, and dt 0.1. 

Figure 1 presents three runs of the AES model. Col- 
umn A displays the dynamics of an hexagonal active crys- 
tal composed of N = 91 agents. At t = (panel Al), 
agents are placed with random orientations on a per- 
fect hexagonal lattice, separated by oIa = 0.65. Nearest 
neighbors are connected by springs with natural length 
I = oIa and spring constant k/l = 5/0.65. Only sensing 
noise is considered here (D r = 0.5, Dq = 0), but results 
remain qualitatively unchanged for other types of noise. 
As time advances, growing regions of coherent motion de- 
velop, eventually deforming the whole structure (A2) un- 
til the group starts translating or rotating collectively. In 
the case shown, the system converges to a rotating state 
where the axis of rotation and barycenter do not coincide 
(A3), thus rotating while translating. Note that rotating 
states will always have higher elastic energy, since inner 
and outer shells cannot move at the same vo speed and 
must be sped up or slowed down by elastic forces. They 
are less frequent and metastable, eventually relaxing to 
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FIG. 1. Active elastic sheet simulations of Eqs. |T]) and |2]). 
A: Hexagonal active crystal at t = (Al), 240 (A2)Jand 1700 
(A3). B: Rod-like active crystal at t = (Bl), 400 (B2), and 
1700 (B3). C: Active solid at same times as column B; darker 
agents symbolize higher local alignment. 



lower-energy, translating solutions. However, we show 
one here to illustrate its dynamics, which cannot be at- 
tained by the Vicsek model. 

Column B displays an active elastic rod, comprised of 
N = 118 agents arranged into three rows, for the same 
noise as in column A. It is generated (Bl) by placing ran- 
domly oriented agents with nearest-neighbor distances 
ds = 0.32 (within rows) and d* B = 0.58 (between rows), 
linking all agents separated by d < 1 with springs of nat- 
ural length d and spring constant k/l = 5/d. Here again, 
larger and larger regions of coherent deformation emerge 
(B2) until CM is attained and the rod starts moving (B3). 
Since the first bending mode has the largest final defor- 
mation, a collective heading direction perpendicular to 
the rod's axis is favored. 

Column C shows TV = 891 agents forming an active 
solid (given the irregular agent positions) with two holes. 
To construct it, we distribute agents at random, homoge- 
neously within the structure, connecting all agents sep- 
arated by d < 1 with springs, following the same pro- 
cedure used in column B. Noise is set here to zero, but 
equivalent dynamics are observed for small enough D r 
and Dq. To highlight ordered regions, each agent's dark- 
ness is displayed proportional to the local order, defined 
similar to i/j but summing only over the focal agent and 
others linked to it, instead of the whole system. Initially, 
most of the structure appears in light gray, since agents 
are randomly oriented (CI). As time advances, growing 
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FIG. 2. (color online). Order parameter and Binder cumu- 
lant vs. positional sensory noise D r and angular actuation 
noise Do for hexagonal active crystals with N — 91 (same as 
on Fig.[l] panel Al) and N — 547 agents. Top panels display 
the mean and local maxima of the distribution of ip values 
obtained in simulations. Insets show these distributions in 
the transition region. For large enough systems, both cases 
display a first order transition with a bistable region. 



regions of coherent motion emerge (C2), until the whole 
structure starts moving when agents become sufficiently 
aligned (C3). 

The AES model displays a discontinuous order- 
disorder transition similar to that in the Vicsek model. 
Figure [2] examines this transition as a function of noise 
for the hexagonal active crystal on column A of Fig.[I]and 
for a larger (N = 547) hexagonal configuration with iden- 
tical parameters. We performed 30 to 80 runs per noise 
value, storing 2000 values of ip per run (every 500 time- 
steps after the initial 10 6 ). Top panels show the mean and 
local maxima of the ip distributions and bottom ones, the 
corresponding Binder cumulants G = 1 — | (^ 4 ) / (^ 2 ) 
[26] . As we increase either sensing noise D r (left) or ac- 
tuation noise Dq (right), ip jumps from an ordered state 
where agents self-organize to a disordered state where 
they continue to point in random directions. The dis- 
played Binder cumulants and bimodal distributions show 
that there is a region of bistability around the critical 
point in both cases. In the sensing noise case, G < 
close to the transition for N = 547, providing evidence 
for bistability in large enough systems. For the actuation 
noise case, N = 547 does not appear to be large enough 
to reach G < 0, but the dip in the transition region drops 
further and further below G = 1/3 (the expected value 
for the disordered phase) as the system size is increased, 
suggesting that G will reach negative values in larger sys- 
tems [I2j [23] [26] . These numerical tests (and others we 
performed with different noise types and agent configu- 
rations) show that the AES transition is first order and 
has a bistable region. 



An interesting aspect of the AES model is that we can 
use a continuous elastic sheet approximation to perform 
analytical calculations. We follow this approach to carry 
out a standard linear stability analysis [24 of the trans- 
lating CM state in the zero noise case. We begin by writ- 
ing the elastic forces F — (F x , F y ) that result from small 
displacements u = (u xi u y ) of points on the membrane 
with respect to their equilibrium positions 



F x = (A • 
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where the elastic constants are given by the Lame pa- 
rameter A and shear modulus \i [24]. We then linearize 
Eqs. ([!]) and Q around an equilibrium solution with un- 
deformed membrane and all agents moving at speed 
in the x direction, obtaining u x = aF Xl ii y = vo<j), and 
(p = pF y , where <\> denotes perturbations to the = 
equilibrium heading angle. Casting these expressions in 
Fourier space with wavevector components (k x ,k y ), we 
can write the perturbation dynamics in matrix form and 
compute its eigenvalues A to determine stability. These 
are found to satisfy the characteristic polynomial equa- 
tion A 3 + C 2 A 2 + CiA + C = 0, with 



C = afinvo (A + 2/i) [k 2 x + k 2 y ] 



d = pvo [fikl + (A - 
C 2 = a [(A + 2/i) k\ 



2/i) k 
-fik 2 ] 
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Using Routh's stability criterion (here given by C1C2 > 
Co [27), we find that the system is stable if 
afivo (A + /i) 2 k 2 k 2 > 0, which is always verified. We con- 
clude that translating CM solutions are always linearly 
stable. This is not the case, however, for most variations 
of the AES model. For example, if we consider a con- 
stant speed algorithm by setting a = 0, the characteristic 
polynomial becomes A 3 + (3 vo [/^k 2 + (A + 2/i) k 2 ] A = 0, 
which only has null or imaginary solutions. Linear per- 
turbations will therefore not dampen out, but produce 
instead permanent oscillations. Numerical simulations 
confirm that, even for zero noise and starting from a per- 
fectly aligned initial condition, the group loses order and 
agents start rotating in place instead of aligning. 

We now characterize the nonlinear energy cascading 
mechanism that drives the AES model to self-organize. 
Figure [3] presents the energy dynamics of an hexagonal 
active crystal simulation with N = 91 and zero noise that 
converges to a translating solution. Top panels display 
the total kinetic and potential energies as a function of 
time. Panel C shows the spectral decomposition of the 
latter into its elastic normal modes, listed in order of 
growing energy and without accounting for degeneracies. 
It is produced by first computing all 182 elastic normal 
modes numerically (without considering agent orienta- 
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FIG. 3. (color online). Kinetic energy (A), elastic energy 
(B), and spectral decomposition of the elastic energy (C) as 
a function of time for an hexagonal N = 91 active crystal 
simulation (with zero noise and same initial condition as in 
panel Al of Fig. [I]) that converges to the aligned, translating 
state. Brighter points on C indicate higher energies. After an 
initial transient, A and B converge to their stationary values 
for collective translational motion. All modes display energy 
levels that oscillate as they decay, with higher modes decaying 
faster. Elastic energy flows to lower modes, producing coher- 
ent motion that eventually reaches the lowest (translational 
or rotational) modes. 



tions or self-propulsion) and then expanding the dynam- 
ics into this basis. The initial condition is set as in panel 
Al of Fig. [T] with zero potential energy and kinetic en- 
ergy Ek = Nvq/2 = 1.82 x 10 -4 (setting the agent mass 
to 1). As the membrane deforms during the initial tran- 
sient, potential energy grows and becomes broadly dis- 
tributed over all modes (as expected for disordered sys- 
tems), while kinetic energy drops. As time advances, the 
system rearranges itself into configurations with lower 
elastic energy and higher kinetic energy, eventually reach- 
ing again (now in the ordered, translating state) values 
close to zero and Ek, respectively. The energy of each 
mode oscillates while decaying, as in an underdamped 
oscillator, with higher modes decaying faster than lower 
ones. This results from a combination of standard elas- 
ticity, self-propulsion, and the coupling between elastic 
forces and turning rate imposed by the AES model. In- 
deed, in standard damped elastic systems, higher energy 
modes also decay faster than lower ones. Here, how- 
ever, each agent is continuously injecting energy through 
its self-propulsion term, so motion cannot be fully damp- 
ened. Instead, modes decay by steering agents away from 



them. Self-propulsion thus feeds energy to lower and 
lower modes, eventually reaching translational or rota- 
tional modes and achieving CM. Note that, despite this 
mechanism, similar models may not converge to CM. 
This will occur, for example, if agents inject too much 
energy into high-energy modes while turning (as in the 
a = case described above) or if they overshoot the an- 
gles that dampen high-energy modes by instantaneously 
switching heading (as in [22j|23]) instead of integrating 
Eq. (§. 

We have identified in this Letter an alternative, 
elasticity-based mechanism that, in contrast to the Vic- 
sek case, requires no exchange of heading information 
to achieve CM. It also requires no switching of interact- 
ing neighbors over time to overcome the Mermin- Wagner 
theorem and achieve long-range order at non-zero noise 
levels [TTJ [I4j [25j [28] . Instead, despite including no fluid- 
like mixing, we found that the AES model displays no loss 
of long-range order for larger structures, although the en- 
ergy cascading mechanism takes longer to converge. We 
found only one other model, introduced in [20] to study 
the collective migration of tissue cells, that can display 
CM under similar conditions. In a version of this algo- 
rithm (designed to study active jamming) where agents 
only have repulsive interactions and are confined to a 
circular box, a similar elasticity-based mechanism was 
shown to be responsible for the dynamics of the jammed 
phase [21]. While that model describes a different situa- 
tion, where interaction forces can displace agents perpen- 
dicular to their heading, its detailed comparison to the 
AES model should yield a more complete understanding 
of the energy cascading mechanism. 

Given that some kind of attraction-repulsion interac- 
tion must be present for any group of moving agents to 
remain cohesive and avoid overlapping, we expect the 
elasticity-based mechanism introduced here to play a rel- 
evant role in the CM dynamics of a variety of systems. 
These could include microscopic biological or robotic 
agents (or even collectively migrating tissue cells [20j [29] ) 
that can exert attraction-repulsion physical forces while 
being too simple to exhibit explicit aligning interac- 
tions. We expect most animal groups to typically com- 
bine elasticity-based and alignment-based mechanisms in 
their CM dynamics, thus effectively behaving as an active 
viscoelastic material composed of aligning spins. An in- 
teresting open question is the extent to which each mech- 
anism is responsible for the CM of specific real-world 
swarms, which can depend on the time-scale and dynam- 
ical state considered. Note that each mechanism results 
from interactions that produce different dynamical signa- 
tures, such as the properties of propagating waves or the 
response to perturbations. Now that there is a growing 
number of experiments that allow the precise tracking of 
different types of swarms [30j [31] , these signatures could 
help determine their individual interaction rules based 
only on observed properties of their collective dynamics. 
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